Emergence of Order in Selection-Mutation Dynamics 
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We characterize the time evolution of a d- dimensional probability distribution by the value of its 
final entropy. If it is near the maximally-possible value we call the evolution mixing, if it is near zero 
we say it is purifying. The evolution is determined by the simplest non-linear equation and contains 
a d x d matrix as input. Since we are not interested in a particular evolution but in the general 
features of evolutions of this type, we take the matrix elements as uniformly-distributed random 
numbers between zero and some specified upper bound. Computer simulations show how the final 
entropies are distributed over this field of random numbers. The result is that the distribution 
crowds at the maximum entropy, if the upper bound is unity. If we restrict the dynamical matrices 
to certain regions in matrix space, for instance to diagonal or triangular matrices, then the entropy 
distribution is maximal near zero, and the dynamics typically becomes purifying. 



I. MOTIVATION 

The abstract setting of evolution equations in fields 
like chemistry, biology, population dynamics, financial 
mathematics and similar subjects can be formulated as 
follows: The state of the system is given by d positive 
numbers pt , which sum to unity and may be interpreted 
as probabilities, as relative populations, or financial as- 
sets. Their time dependence is governed by a dynamical 
equation, which expresses the time change of pi by the 
value of all the pj . The basic question is whether in the 
course of time the state will purify (which means that one 
Pi will tend to unity and all the others to zero) or will 
become completely mixed (all pi become equal to 1 /d) . 
To measure this tendency, we use the quadratic entropy 



s = ^2pi(i -Pi)- 



(i) 



S vanishes for a pure state, it is (d — l)/d for the com- 
pletely mixed state, and is in between otherwise. The 
evolution equations are constructed such that the ba- 
sic properties of the pi are preserved and their values 
asymptotically tend to a limit independent of the initial 
conditions. However, this limit will depend on the coeffi- 
cients oiij appearing in the evolution equation as shown 
below. They are supposed to be determined by the ac- 
cidental circumstances of the system and should be con- 
sidered as random numbers. The basic issue therefore 
is, whether the asymptotic values for S, taken as a func- 
tion of these coefficients, tend to cluster around zero, or 
around (d — l)/d. Correspondingly, we say that these 
evolution equations are either purifying or mixing. 

We answer this question for the simplest type of evo- 
lution equation which, speaking financially, says the fol- 
lowing: per unit time each player i receives the fraction 



o-ij of the assets of the players j. Then, there is an overall 
tax such that the total amount of money remains unity. 
If we take all the ay to be random numbers between 
and 1, we find that the asymptotic values for S cluster 
just below (d — and this distribution is the sharper 
the bigger the dimension d. Only for special classes of 
ctij they cluster near zero. This means that one can 
make a fortune only by an intelligent strategy, or that 
the "survival of the fittest" only holds under special cir- 
cumstances. 

The paper is organized as follows: In Section |TT] wc 
define our model and outline the numerical procedure 
for the computation of the entropy S. In Section IIIII wc 
summarize our numerical results for three to eight dimen- 
sional systems. The qualitative behavior is already ap- 
parent for the two-dimensional case, which may be solved 
analytically. This is demonstrated in Sect ion HVl We con- 
clude in Section [V] with some remarks on fitness and the 
quantum- mechanical generalization. 



II. DYNAMICAL EVOLUTION INCLUDING 
MUTATION 

We consider a d-dimensional dynamical system with a 
state vector 

p = {pi}, < pi < 1; i = 1, 2, . . . , d, 

which evolves according to 



, d d d 

— = a ljPj - pt 2^ 2^ a i kPk 
j=l j=i fe=i 



(2) 



The coefficients cty are elements of a d x d matrix a and 
are assumed to be strictly positive random variables with 
upper bounds specified below. From the general solution 
ofEq. @, 
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Pit) 



exp (at)p(0) 
EiLjexp (at)p(0)]i 



(3) 
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it follows that the vector p is constrained to the simplex 



= 1 . 



(4) 



i=i 



provided the initial point is also contained in Sd, 
£i=iPi(P) = 1- The state p may be regarded as a d- 
dimensional probability distribution. The dynamics of 
Equation ([2]) is closely related to the quasi-species equa- 
tion introduced by Eigen and Schuster [2j , which has been 
extensively studied in the past [|[ . It also corresponds to 
the evolution equation introduced by Nowak, Komarova 
and Niyogi for the study of language learning by children 
In terms of chemical reaction kinetics, it models evo- 
lutionary dynamics based on selection and mutation. 

a is considered a matrix and, therefore, is similar to 
an orthogonal sum of matrices in a Jordan normal form, 

a = M- 1 J2( a i + S i) M ^ 



where the sum is over all Jordan blocks corresponding to 
different eigenvalues of a with degeneracy dj, and S is 
the shift, 



S = 



/0 1 
10 
1 

V . . . 



(•5) 



The matrix a is diagonal in each dj x dj Jordan block. 
exp(a + S) can easily be calculated, because a and S 
commute and S d — 0. Since the off-diagonal elements of 
a are assumed here to be strictly positive, the Perron- 
Frobenius theorem applies. As a consequence, the max- 
imum eigenvalue, A, of a is always non degenerate and 
strictly positive, A > 0. The corresponding eigenvec- 
tor, v, is real and positive. The associated Jordan block 
is one-dimensional. It follows that the leading term of 
the numerator of Eq. ([3]), for t — > oo, is e xt v. The 
exponential is cancelled by an analogous expression for 
the denominator. We conclude from Eq. ^ that the 
asymptotic equilibrium is determined by the eigenvector 
v belonging to A and is normalized according to Eq. ^ 
and, thus, confined to the simplex Sd' 



lim p{t) 




(6) 



This asymptotic solution is simply denoted by p for the 
remainder of this paper. 

To quantify the asymptotic behavior of the system, we 
use the quadratic entropy defined in Eq. ([1]), which is 
bounded, 



0<S< 



1 



S 150 




FIG. 1: (Color online) Random matrices with D = U — L — 
1: Normalized entropy distributions p(S) for various dimen- 
sions d as specified by the labels. 



The upper bound corresponds to a maximally mixed 
state, where all the pi are equal, pi = l/d. We say the 
dynamics is maximally mixing, if this state is approached 
for t — > oo. The lower bound of S corresponds to a pure 
state, where only one component of p approaches unity 
and all the others vanish. Accordingly, the dynamics is 
said to be maximally purifying, if such a state is asymp- 
totically reached. We say the system is (predominantly) 
mixing or purifying, if the distribution p(S) is peaked at 
or near the respective upper or lower bound. 



III. RESULTS 

In the following we consider classes of random matri- 
ces a with strictly positive elements, and determine the 
asymptotic behavior of the system by computing the nor- 
malized entropy distribution p(S), ji d p(S)dS = 1, 
for each class. This is achieved by numerically diag- 
onalyzing (typically 10 7 to 10 s ) reaction-rate matrices 
a and constructing the entropy histogram from the re- 
normalized eigenvectors, Eq. belonging to the re- 
spective maximum eigenvalues. For the classification of 
a we distinguish between i) diagonal elements (D), ii) el- 
ements above the diagonal (U), and iii) elements below 
the diagonal (L). The elements are randomly drawn from 
a uniform distribution such that an upper bound exists 
for each group: 



< a u < D ; 1 < i < d 

< OLij < U ; 1 < i < j < d 

< Q-ij < L ; 1 < j < i < d 



(8) 
(9) 



(7) 

We specify various classes of a by D, U, and L. 
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FIG. 2: (Color online) Entropy distributions p(S) in three 
dimensions (d — 3) for random matrices with D — 1 and 
upper bounds for the off-diagonal elements U — L varying 
between 10 -4 and unity. The distributions are clipped at 40 
for display purposes. The selected off-diagonal bounds U = L 
are arranged on a logarithmic scale. 



FIG. 3: (Color online) Entropy distributions p(S) in eight 
dimensions (d — 8) for random matrices with D = 1 and 
upper bounds for the off-diagonal elements U — L varying 
between 10 -4 and unity. The distributions are clipped for 
p > 40. The selected off-diagonal bounds U = L are arranged 
on a logarithmic scale. 



A. The most random case D — U = L = 1 

First we consider in Fig. [1] the most random case, 
where the upper bounds for all matrix elements are the 
same, D = U = L = 1. The dimension d is varied 
from 2 to 10. Not surprisingly, the off-diagonal elements 
are responsible for effective mixing, and the distributions 
p(S) are sharply peaked near their most-random edge, 
(d— l)/d. For d > 3 they become narrower for increasing 
d, and are expected to approach <5(1) for d — > oo. 



B. The case D = 1 and U = L > 

In Fig. [2]the diagonal elements for a three-dimensional 
system are bounded by unity, D = 1. Various entropy 
distributions are shown, for which the bounds for the 
off-diagonal elements are equal, U = L, and are varied 
from 1CP 4 to 10 4 . Note the logarithmic scale. The 
system is clearly purifying if the off-digonal elements are 
small. However, it becomes mxing for U = L > 0.1. This 
transition takes place for smaller values of U = L, if the 
dimension of the system is larger. This is demonstrated 
in Fig. OS for d=8. 



C. Strong asymmetry: The dependence on D and 
L for U — 1 

We stick with three dimensions and consider the case 
where the bound for the upper off-diagonal elements is 
fixed, U = 1, whereas the respective bound for the lower 
off-diagonal elements is varied in steps over a huge range, 
10~ 5 < L < 10 5 . The top and bottom figures of Fig. g] 
differ by the choice of D, the upper bound of the diagonal 



elements. In the top figure, D is set to vanish, whereas in 
the figure at the bottom D is set to unity. Although the 
off-diagonal elements may be small, they are still strictly 
positive and the conditions of Section |TT] still apply. One 
infers from the top of Fig. 2] that for vanishing D the 
system - very gradually - becomes more purifying, if the 
ratio of upper bounds U/L becomes large or small. In the 
limit L = the system is strictly purifying, although the 
convergence towards this state is slow as we have verified 
by direct integration of the equations of motion @ . 

The entropy distributions very often have peaks which 
belong to mixed substructures with a lower dimension 
(such as at S = 0.5 in Fig. rj|. This feature may be 
observed with many other examples given in this paper. 

To complement these results for three-dimensional sys- 
tems, we show in Figure 03 how the shape of the entropy 
distribution depends on the upper bound of the diagonal 
elements D. Again, U = 1, whereas the upper bound 
for the lower off-diagonal elements, L, is 10~ 2 , 10~ 4 , and 
10~ 6 for the figures from top to bottom, respectively. 
The surprising result is that for already rather small (but 
strictly positive) matrix elements below the diagonal ( 
for L = 10~ 2 at the top of Fig. 03), the system is mix- 
ing for small and vanishing (not shown) D, where the 
rather broad maximum of the entropy distribution only 
very slowly moves towards the purifying case S = for 
much smaller values of L. (middle and bottom figure). 

This results also applies for higher dimensions. In Fig. 
03 we show the eight-dimensional analogue to Fig. 03 
where again p(S) is shown for various D, where U is 
unity, and L varies from 10~ 2 to 10~ 8 for the figures 
from top to bottom. Only for U — the system becomes 
strictly purifying for small D. 
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d = 3: D = Q,U= 1 




FIG. 4: (Color online) Entropy distributions p(S) in three 
dimensions (d — 3) for various L (logarithmic scale). The 
bound U = 1, whereas D for the diagonal elements vanishes 
at the top, and is unity at the bottom. The distributions are 
clipped for p > 40. 



IV. THE TWO-DIMENSIONAL CASE 

In an attempt to understand some of the details of this 
complicated behavior, we turn to the two-dimensional 
case, d = 2, which may be solved analytically, at least in 
principle. The behavior in two dimensions is qualitatively 
the same as for higher dimensions. This is demonstrated 
by the simulation results in Fig. where entropy den- 
sities are shown for various D varying over many orders 
of magnitude, and for bounds L, decreasing, from top to 
bottom, from 10 -3 to 1CP 5 , respectively. S is bounded 
between and 0.5 in this case. 

The 2x2 matrix of rate constants is written 
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FIG. 5: (Color online) Entropy distributions p(S) in three 
dimensions (d = 3) for various D (logarithmic scale). The 
bound for the upper off-diagonal elements U = 1, whereas 
the respective bounds for the lower off-diagonal elements are 
L — 10~ 2 , 10~ 4 , and 10~ 6 for the top, middle and bottom 
figures, respectively. The distributions are clipped for p > 40. 



The maximum eigenvalue of a is 



a b 
c d 



(10) 



where, in accordance with Eq. [9j the strictly positive 
elements are taken to be bounded by 



< {a,d} < D ; <6 <U ; <c < L. 



(11) 




According to Sec. UH the corresponding (properly nor- 
malized) eigenvector p with components (pijPa) gives the 
asymptotic state of the system. From the eigenvalue 
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FIG. 6: (Color online) Entropy distributions p(S) in eight 
dimensions (d = 8) for various D (logarithmic scale). The 
bound for the upper off-diagonal elements U = 1, whereas 
the respective bounds for the lower off-diagonal elements are, 
from top to bottom, L = 1CT 2 , 10" 4 , 10" 6 , and 1CT 8 . The 
distributions are clipped for p > 40. 



FIG. 7: (Color online) Entropy distributions p(S) in two di- 
mensions (d — 2) for various D (logarithmic scale). The 
bound for the upper off-diagonal element U = 1, whereas 
the respective bounds for the lower off-diagonal element are 
L — 10~ 2 , 10~ 3 , and 1CP 4 for the top, middle and bottom 
figures, respectively. The distributions are clipped for p > 40. 



equation ap = Xp it immediately follows that the ratio 



Pi 

P2 



V/? 2 + be + /3 



(12) 



with f3 = (a — d)/2, only depends on the difference of the 
diagonal elements. The bounds for (3 follow from those 
for {a, 6}, Eq. [□] -D/2 < (3 < D/2. Since a and d are 
assumed to be uniformly distributed, (3 has a triangular 



6 



distribution, 



If 2 \ f D ' 2 



(P)df3 = 1. (13) 



The ratio r — px/p2 is another indicator for mixing: 
A pure state corresponds to r — and r = oo, and the 
maximally mixed state has r = 1. r is related to the 
quadratic entropy by 



S(r) 



2r 



(1 + r) 2 



s 



Once the normalized distribution for r, p(r), is known, 
the corresponding entropy distribution p(S) follows from 



g 
2 




p(S) 



gHg)) 

|dS(r)/dr| 



2S*V1 - 25 



p(r+(5)) 



p(r_(S)).(14) 



Here, r+ is used for the branch 1 < r < oo, and r_ for 
the branch < r < 1: 

1 



r±(S) = - (1- S±Vl-2S 



The distribution for r is obtained from 
p(r) 



1 



(7 /-L /.D/2 

d& / cfc / d(3 

JO J-D/2 



xtt(/3)S r 



(15) 



Here, the <5-function takes care of the constraining rela- 
tion Eq. [12] between r and the matrix elements b, c and 
P, and may be rewritten according to 



Sir- 



V/? 2 + be + \ _ 1_ 

c / 2r V r 



Without loss of generality we consider r = r + > 1 Since 

' dpf(0)*(P -x) = f(x)G - |. 

_D/2 V Z 

where Q is the unit step function, the distribution for r 
becomes 



p(r) 



1 



db / <ic 



b 

rc 

r 



With the transformation 

x = b/Dr ; y = cr/D 



FIG. 8: (Color online) Entropy distributions p(S) in two di- 
mensions (d = 2) for U = 1 and D = 10" 3 . As indicated 
by the labels, the two distributions are for L = 10 -4 (also 
marked with an arrow in the bottom figure of Fig. [7} and 
L = 1CT 5 . The points are simulation results, and the smooth 
lines are the theoretical predictions as described in the main 
text. 



one obtains 

p(r) = {rxayay 1 I dxl dy(l - \y - x\){x + y) 



/0 ^0 

xe(l-\y-x\), 



(16) 



where xq = U/Dr and yo = Lr/D. The domain for 
the remaining integration in the xy plane is the intersec- 
tion of a rectangle, defined by < x < xo; < y < yo, 
and of a strip given by— 1 + x < y < 1 + x. For the 
general case this domain is slightly involved. Here, we 
confine ourselves to a specific situation, which is readily 
solved. For example, we take U = 1, L = 10~ 4 , and 
D = 1CP 3 . The simulation results for this particular 
case are marked by an arrow in the bottom panel of Fig. 
[7] Under these conditions, p(r) becomes a maximum for 
r = ?*o = y/U/L > 1, for which xo — yo and the overlap 
of the rectangle with the strip is maximized. 

For the regime of very large r, r = r+ > ro > 1, a di- 
rect evaluation of Eq. (fl~6]) gives p(r) = (U/L)r~ 3 . Since 
only r+ contributes, the entropy distribution is given by 
the first term of Eq. (fT4|) . To lowest order in S, one 
finds p(S) w (U/AL)S for S < S = 2r /(l +r ) 2 , which 
agrees well with the numerical results. This is shown in 
Fig.! 

For the regime of slightly smaller r, ro > r = r+ > 1, 
a similar evaluation of the integrals in Eq. (|16|) yields 
p(r) = {L/U)r and, hence, p(S) w (<iL/U)/S 3 for 
S > Sq. Again, this is in complete agreement with the 
numerical data (Fig. [8|). 

We conclude that all the details of the two-dimensional 
entropy distributions in Figure [7] are accessible to ana- 
lytical computation. This is of some significance because 
of the qualitative similarity with the results for higher 
dimensions. 
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V. FITNESS 
A. The classical case 
If the evolution equations of Eq. [2] are written as 
Pi =^2 a ijP] ~ #i> 

3 

the global fitness of the state, 

j k 

evolves according to 

<j> = ^ ^ (Xjk Ct M Pl ~ Pk ^ ^ CtlmPm 
j k \ I m / 

It increases if a is diagonal. 

B. The quantum setting 

Our results so far can readily be generalized to the 
quantum mechanical setting [8j, where the distribution 
Pl,P2,••• ,Pd is replaced by a density matrix p which 
obeys 

dp(t)/dt = hp{t) + p(t)h 1 - p{t)Tr(hp(t) + p(i)rf). (17) 

h is a sort of Hamilton operator, and is its Hermi- 
tian adjoint. The solution is almost the same as in the 
classical case before, 

p(t) = exp(M)p(0)exp(i/i t )[Trexp(M)p(0)exp(</i t )]" 1 , 

just observe the non-commutativity of the operators. En- 
tropy S and fitness F generalize to 

S = Trp(l - p); F = Trp(h + = dc/dt, 

with 

c(t) = ln[Trexp(M)p(0)exp(i/i t )]. 

Thus, the increase in fitness means convexity of the func- 
tion c(t). The increase of F for diagonal a in the classical 
case becomes now the following 
Proposition: 

F is increasing if h is Hermitian. 
Proof: 

The evolution equation (117|) tells us generally 

dF/dt = Trp{h+tf) = Tr[hp+prf-pTrp(h+rf)}{h+tf). 

For Hermitian h this becomes four times Tr(h 2 p) — 
(Trph) 2 . The generalization of the Schwartz inequality 
to traces of operator products implies 

(Trph) 2 = (Tr^/y/ 2 ^) 2 < (Trp)(Trph 2 ), 



which gives F < for all p. 
Remarks: 

1) Without Hermiticity, F need not increase. This ap- 
pears already for a nilpotent h, h = 0. Then exp(fti) = 
1 + hi and 

ln[Tr(l + ht)p{\ + th r )} = lnTr[(l + th + th) + t 2 hh))p] 

is no longer convex in t. 

2) The Schwarz inequality becomes an equality only 
if the two vectors are proportional. Therefore the proof 
tells us that the fitness keeps increasing unless p projects 
onto an eigenvector of h. Of these, only the one with the 
highest eigenvalue has the biggest F and, therefore, is an 
attractor. Excluding degeneracy this means that then 
also the entropy goes to zero and the system becomes 
pure. S cannot be a Lyapunov function, because near a 
repellor it will be very small, so it has to increase only to 
decrease again near the attractor. 

3) The situation is just the converse for the Lindblad 
dynamics [5j], which provides a consistent description of 
a dissipative quantum system [1, 0, Q : 

p = hph) - - (hrfp + phtf) . (18) 

This equation is linear in p but quadratic in h. There, 
the fitness F = Trhph^ increases for a nilpotent h, F — 
h^hphh^ > 0, and is constant for h — . 

VI. CONCLUDING REMARKS 

The motivation for this study was to see, whether 
the nonlinear evolution equations let order emerge from 
chaos purely by accident. The answer depends where we 
let chance act. Our results can be simply stated as fol- 
lows: According to Eq. ((9|) the matrix space is spanned 
by three parts, the diagonal matrices (D), and the upper 
(U) and lower (L) triangular matrices. Matrices taken 
from solely one of these subalgebras give a purifying dy- 
namics, but admixture of a small fraction of another part 
renders the dynamics mixing. Since zero is not a random 
number, these subalgebras are never pure, but we keep 
the admixtures so small that they do not matter. These 
features can be easily understood in the financial setting. 
There, the pi are the assets of party i, and the rules of 
the game are that per unit time party k gives the ratio 
otik of its asset to party i. This money is not deducted 
from its account since the aik are strictly positive, but, 
instead, there is a flat tax which keeps the total amount 
of money fixed. The diagonal elements correspond to the 
case where each party makes its own investments and 
eventually the one with the highest interest rate becomes 
the richest. The triangular dynamics means that there 
is an order between the parties, and money is handed 
only from parties of lower order to the ones of higher 
order. Obviously, the one on top of this order pyramid 
eventually wins. 
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In a mixed dynamics these tendencies will compete, 
and since there are so many more possibilities for them 
to act in different directions, a random choice of mixtures 
will lead to a mixing dynamics. 

Thus, within this equation only the simplest and most 
brutal strategies lead to the dominance of one party, but 
fortuna may easily even this out. 
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